home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / C / Applications / POV-Ray 3.0.2 / src / SOURCE / PRISM.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-10-09  |  39.2 KB  |  1,897 lines  |  [TEXT/CWIE]
  1. /****************************************************************************
  2. *                   prism.c
  3. *
  4. *  This module implements functions that manipulate prisms.
  5. *
  6. *  This module was written by Dieter Bayer [DB].
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file. If
  16. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  18. *  Forum.  The latest version of POV-Ray may be found there as well.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. /****************************************************************************
  27. *
  28. *  Explanation:
  29. *
  30. *    The prism primitive is defined by a set of points in 2d space which
  31. *    are interpolated by linear, quadratic, or cubic splines. The resulting
  32. *    2d curve is swept along a straight line to form the final prism object.
  33. *
  34. *    All calculations are done in the object's (u,v,w)-coordinate system
  35. *    with the (w)-axis being the sweep axis.
  36. *
  37. *    One spline segment in the (u,v)-plane is given by the equations
  38. *
  39. *      fu(t) = Au * t * t * t + Bu * t * t + Cu * t + Du  and
  40. *      fv(t) = Av * t * t * t + Bv * t * t + Cv * t + Dv,
  41. *
  42. *    with the parameter t ranging from 0 to 1.
  43. *
  44. *    To intersect a ray R = P + k * D transformed into the object's
  45. *    coordinate system with the linear swept prism object, the equation
  46. *
  47. *      Dv * fu(t) - Du * fv(t) - (Pu * Dv - Pv * Du) = 0
  48. *
  49. *    has to be solved for t. For valid intersections (0 <= t <= 1)
  50. *    the corresponding k can be calculated from equation
  51. *
  52. *      Pu + k * Du = fu(t) or Pv + k * Dv = fv(t).
  53. *
  54. *    In the case of conic sweeping the equation
  55. *
  56. *      (Pv * Dw - Pw * Dv) * fu(t) - (Pu * Dw - Pw * Du) * fv(t)
  57. *                                  + (Pu * Dv - Pv * Du) = 0
  58. *
  59. *    has to be solved for t. For valid intersections (0 <= t <= 1)
  60. *    the corresponding k can be calculated from equation
  61. *
  62. *      Pu + k * Du = (Pw + k * Dw) * fu(t) or
  63. *      Pv + k * Dv = (Pw + k * Dw) * fv(t).
  64. *
  65. *    Note that the polynomal to solve has the same degree as the
  66. *    spline segments used.
  67. *
  68. *    Note that Pu, Pv, Pw and Du, Dv, Dw denote the coordinates
  69. *    of the vectors P and D.
  70. *
  71. *  Syntax:
  72. *
  73. *    prism {
  74. *      [ linear_sweep | cubic_sweep ]
  75. *      [ linear_spline | quadratic_spline | cubic_spline ]
  76. *
  77. *      height1, height2,
  78. *      number_of_points
  79. *
  80. *      <P(0)>, <P(1)>, ..., <P(n-1)>
  81. *
  82. *      [ open ]
  83. *      [ sturm ]
  84. *    }
  85. *
  86. *    Note that the P(i) are 2d vectors.
  87. *
  88. *  ---
  89. *
  90. *  Ideas for the prism was taken from:
  91. *
  92. *    James T. Kajiya, "New techniques for ray tracing procedurally
  93. *    defined objects", Computer Graphics, 17(3), July 1983, pp. 91-102
  94. *
  95. *    Kirk ...
  96. *
  97. *  ---
  98. *
  99. *  May 1994 : Creation. [DB]
  100. *
  101. *****************************************************************************/
  102.  
  103. #include "frame.h"
  104. #include "povray.h"
  105. #include "vector.h"
  106. #include "povproto.h"
  107. #include "bbox.h"
  108. #include "matrices.h"
  109. #include "objects.h"
  110. #include "polysolv.h"
  111. #include "prism.h"
  112.  
  113.  
  114.  
  115. /*****************************************************************************
  116. * Local preprocessor defines
  117. ******************************************************************************/
  118.  
  119. /* Minimal intersection depth for a valid intersection. */
  120.  
  121. #define DEPTH_TOLERANCE 1.0e-4
  122.  
  123. /* |Coefficients| < COEFF_LIMIT are regarded to be 0. */
  124.  
  125. /*#define COEFF_LIMIT 1.0e-20 */
  126.  
  127. #define COEFF_LIMIT 1.0e-16 /*changed by CEY 11/18/95 */
  128.  
  129. /* Part of the prism hit. */
  130.  
  131. #define BASE_HIT   1
  132. #define CAP_HIT    2
  133. #define SPLINE_HIT 3
  134.  
  135.  
  136.  
  137. /*****************************************************************************
  138. * Local typedefs
  139. ******************************************************************************/
  140.  
  141.  
  142.  
  143. /*****************************************************************************
  144. * Static functions
  145. ******************************************************************************/
  146.  
  147. static int intersect_prism PARAMS((RAY *Ray, PRISM *Prism, PRISM_INT *Intersection));
  148. static int in_curve PARAMS((PRISM *Prism, DBL u, DBL v));
  149. static int test_rectangle PARAMS((VECTOR P, VECTOR D, DBL x1, DBL y1, DBL x2, DBL y2));
  150. static int   All_Prism_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  151. static int   Inside_Prism PARAMS((VECTOR point, OBJECT *Object));
  152. static void  Prism_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  153. static void  *Copy_Prism PARAMS((OBJECT *Object));
  154. static void  Translate_Prism PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  155. static void  Rotate_Prism PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  156. static void  Scale_Prism PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  157. static void  Transform_Prism PARAMS((OBJECT *Object, TRANSFORM *Trans));
  158. static void  Invert_Prism PARAMS((OBJECT *Object));
  159. static void  Destroy_Prism PARAMS((OBJECT *Object));
  160.  
  161.  
  162. /*****************************************************************************
  163. * Local variables
  164. ******************************************************************************/
  165.  
  166. static METHODS Prism_Methods =
  167. {
  168.   All_Prism_Intersections,
  169.   Inside_Prism, Prism_Normal,
  170.   Copy_Prism,
  171.   Translate_Prism, Rotate_Prism,
  172.   Scale_Prism, Transform_Prism, Invert_Prism, Destroy_Prism
  173. };
  174.  
  175.  
  176.  
  177. /*****************************************************************************
  178. *
  179. * FUNCTION
  180. *
  181. *   All_Prism_Intersections
  182. *
  183. * INPUT
  184. *
  185. *   Object      - Object
  186. *   Ray         - Ray
  187. *   Depth_Stack - Intersection stack
  188. *   
  189. * OUTPUT
  190. *
  191. *   Depth_Stack
  192. *   
  193. * RETURNS
  194. *
  195. *   int - TRUE, if a intersection was found
  196. *   
  197. * AUTHOR
  198. *
  199. *   Dieter Bayer
  200. *   
  201. * DESCRIPTION
  202. *
  203. *   Determine ray/prism intersection and clip intersection found.
  204. *
  205. * CHANGES
  206. *
  207. *   May 1994 : Creation.
  208. *
  209. ******************************************************************************/
  210.  
  211. static int All_Prism_Intersections(Object, Ray, Depth_Stack)
  212. OBJECT *Object;
  213. RAY *Ray;
  214. ISTACK *Depth_Stack;
  215. {
  216.   int i, max_i, Found;
  217.   PRISM_INT *Inter;
  218.   VECTOR IPoint;
  219.  
  220.   Found = FALSE;
  221.  
  222.   Inter = ((PRISM *)Object)->Intersections;
  223.  
  224.   max_i = intersect_prism(Ray, (PRISM *)Object, Inter);
  225.  
  226.   if (max_i)
  227.   {
  228.     for (i = 0; i < max_i; i++)
  229.     {
  230.       if ((Inter[i].d > DEPTH_TOLERANCE) && (Inter[i].d < Max_Distance))
  231.       {
  232.         VEvaluateRay(IPoint, Ray->Initial, Inter[i].d, Ray->Direction);
  233.  
  234.         if (Point_In_Clip(IPoint, Object->Clip))
  235.         {
  236.           push_entry_i1_i2_d1(Inter[i].d,IPoint,Object,Inter[i].t,Inter[i].n,Inter[i].w,Depth_Stack);
  237.  
  238.           Found = TRUE;
  239.         }
  240.       }
  241.     }
  242.   }
  243.  
  244.   return(Found);
  245. }
  246.  
  247.  
  248.  
  249. /*****************************************************************************
  250. *
  251. * FUNCTION
  252. *
  253. *   intersect_prism
  254. *
  255. * INPUT
  256. *
  257. *   Ray   - Ray
  258. *   Prism - Prism
  259. *   Int   - Prism intersection structure
  260. *   
  261. * OUTPUT
  262. *
  263. *   Int
  264. *   
  265. * RETURNS
  266. *
  267. *   int - Number of intersections found
  268. *   
  269. * AUTHOR
  270. *
  271. *   Dieter Bayer
  272. *   
  273. * DESCRIPTION
  274. *
  275. *   Determine ray/prism intersection.
  276. *
  277. *   NOTE: Order reduction cannot be used.
  278. *
  279. * CHANGES
  280. *
  281. *   May 1994 : Creation.
  282. *
  283. ******************************************************************************/
  284.  
  285. static int intersect_prism(Ray, Prism, Intersection)
  286. RAY *Ray;
  287. PRISM *Prism;
  288. PRISM_INT *Intersection;
  289. {
  290.   int i, j, n;
  291.   DBL k, u, v, w, h, len;
  292.   DBL x[4], y[3];
  293.   DBL k1, k2, k3;
  294.   VECTOR P, D;
  295.   PRISM_SPLINE_ENTRY Entry;
  296.  
  297.   /* Don't test degenerate prisms. */
  298.  
  299.   if (Test_Flag(Prism, DEGENERATE_FLAG))
  300.   {
  301.     return(FALSE);
  302.   }
  303.  
  304.   Increase_Counter(stats[Ray_Prism_Tests]);
  305.  
  306.   /* Init intersection structure. */
  307.  
  308.   Intersection->d =
  309.   Intersection->w = 0.0;
  310.   Intersection->n =
  311.   Intersection->t = 0;
  312.  
  313.   /* Transform the ray into the prism space */
  314.  
  315.   MInvTransPoint(P, Ray->Initial, Prism->Trans);
  316.  
  317.   MInvTransDirection(D, Ray->Direction, Prism->Trans);
  318.  
  319.   VLength(len, D);
  320.  
  321.   VInverseScaleEq(D, len);
  322.  
  323.   /* Test overall bounding rectangle. */
  324.  
  325. #ifdef PRISM_EXTRA_STATS
  326.   Increase_Counter(stats[Prism_Bound_Tests]);
  327. #endif
  328.  
  329.   if (((D[X] >= 0.0) && (P[X] > Prism->x2)) ||
  330.       ((D[X] <= 0.0) && (P[X] < Prism->x1)) ||
  331.       ((D[Z] >= 0.0) && (P[Z] > Prism->y2)) ||
  332.       ((D[Z] <= 0.0) && (P[Z] < Prism->y1)))
  333.   {
  334.     return(FALSE);
  335.   }
  336.  
  337. #ifdef PRISM_EXTRA_STATS
  338.   Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  339. #endif
  340.  
  341.   /* Number of intersections found. */
  342.  
  343.   i = 0;
  344.  
  345.   /* What kind of sweep is used? */
  346.  
  347.   switch (Prism->Sweep_Type)
  348.   {
  349.     /*************************************************************************
  350.     * Linear sweep.
  351.     **************************************************************************/
  352.  
  353.     case LINEAR_SWEEP :
  354.  
  355.       if (fabs(D[Y]) < EPSILON)
  356.       {
  357.         if ((P[Y] < Prism->Height1) || (P[Y] > Prism->Height2))
  358.         {
  359.           return(FALSE);
  360.         }
  361.       }
  362.       else
  363.       {
  364.         if (Test_Flag(Prism, CLOSED_FLAG))
  365.         {
  366.           /* Intersect ray with the cap-plane. */
  367.  
  368.           k = (Prism->Height2 - P[Y]) / D[Y];
  369.  
  370.           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  371.           {
  372.             u = P[X] + k * D[X];
  373.             v = P[Z] + k * D[Z];
  374.  
  375.             if (in_curve(Prism, u, v))
  376.             {
  377.               Intersection[i].t   = CAP_HIT;
  378.               Intersection[i++].d = k / len;
  379.             }
  380.           }
  381.  
  382.           /* Intersect ray with the base-plane. */
  383.  
  384.           k = (Prism->Height1 - P[Y]) / D[Y];
  385.  
  386.           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  387.           {
  388.             u = P[X] + k * D[X];
  389.             v = P[Z] + k * D[Z];
  390.  
  391.             if (in_curve(Prism, u, v))
  392.             {
  393.               Intersection[i].t   = BASE_HIT;
  394.               Intersection[i++].d = k / len;
  395.             }
  396.           }
  397.         }
  398.       }
  399.  
  400.       /* Intersect ray with all spline segments. */
  401.  
  402.       if ((fabs(D[X]) > EPSILON) || (fabs(D[Z]) > EPSILON))
  403.       {
  404.         for (j = 0; j < Prism->Number; j++)
  405.         {
  406.           Entry = Prism->Spline->Entry[j];
  407.  
  408.           /* Test spline's bounding rectangle (modified Cohen-Sutherland). */
  409.  
  410. #ifdef PRISM_EXTRA_STATS
  411.           Increase_Counter(stats[Prism_Bound_Tests]);
  412. #endif
  413.  
  414.           if (((D[X] >= 0.0) && (P[X] > Entry.x2)) ||
  415.               ((D[X] <= 0.0) && (P[X] < Entry.x1)) ||
  416.               ((D[Z] >= 0.0) && (P[Z] > Entry.y2)) ||
  417.               ((D[Z] <= 0.0) && (P[Z] < Entry.y1)))
  418.           {
  419.             continue;
  420.           }
  421.  
  422.           /* Number of roots found. */
  423.  
  424.           n = 0;
  425.  
  426.           switch (Prism->Spline_Type)
  427.           {
  428.             case LINEAR_SPLINE :
  429.  
  430. #ifdef PRISM_EXTRA_STATS
  431.               Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  432. #endif
  433.  
  434.               /* Solve linear equation. */
  435.  
  436.               x[0] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
  437.  
  438.               x[1] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
  439.  
  440.               if (fabs(x[0]) > EPSILON)
  441.               {
  442.                 y[n++] = -x[1] / x[0];
  443.               }
  444.  
  445.               break;
  446.  
  447.             case QUADRATIC_SPLINE :
  448.  
  449. #ifdef PRISM_EXTRA_STATS
  450.               Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  451. #endif
  452.  
  453.               /* Solve quadratic equation. */
  454.  
  455.               x[0] = Entry.B[X] * D[Z] - Entry.B[Y] * D[X];
  456.  
  457.               x[1] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
  458.  
  459.               x[2] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
  460.  
  461.               n = Solve_Polynomial(2, x, y, FALSE, 0.0);
  462.  
  463.               break;
  464.  
  465.             case CUBIC_SPLINE :
  466.  
  467.               if (test_rectangle(P, D, Entry.x1, Entry.y1, Entry.x2, Entry.y2))
  468.               {
  469. #ifdef PRISM_EXTRA_STATS
  470.                 Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
  471. #endif
  472.  
  473.                 /* Solve cubic equation. */
  474.  
  475.                 x[0] = Entry.A[X] * D[Z] - Entry.A[Y] * D[X];
  476.  
  477.                 x[1] = Entry.B[X] * D[Z] - Entry.B[Y] * D[X];
  478.  
  479.                 x[2] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
  480.  
  481.                 x[3] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
  482.  
  483.                 n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
  484.               }
  485.  
  486.               break;
  487.           }
  488.  
  489.           /* Test roots for valid intersections. */
  490.  
  491.           while (n--)
  492.           {
  493.             w = y[n];
  494.  
  495.             if ((w >= 0.0) && (w <= 1.0))
  496.             {
  497.               if (fabs(D[X]) > EPSILON)
  498.               {
  499.                 k = (w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X] - P[X]) / D[X];
  500.               }
  501.               else
  502.               {
  503.                 k = (w * (w * (w * Entry.A[Y] + Entry.B[Y]) + Entry.C[Y]) + Entry.D[Y] - P[Z]) / D[Z];
  504.               }
  505.  
  506.               /* Verify that intersection height is valid. */
  507.  
  508.               h = P[Y] + k * D[Y];
  509.  
  510.               if ((h >= Prism->Height1) && (h <= Prism->Height2))
  511.               {
  512.                 Intersection[i].t   = SPLINE_HIT;
  513.                 Intersection[i].n   = j;
  514.                 Intersection[i].w   = w;
  515.                 Intersection[i++].d = k / len;
  516.               }
  517.             }
  518.           }
  519.         }
  520.       }
  521.  
  522.       break;
  523.  
  524.     /*************************************************************************
  525.     * Conic sweep.
  526.     **************************************************************************/
  527.  
  528.     case CONIC_SWEEP :
  529.  
  530.       if (fabs(D[Y]) < EPSILON)
  531.       {
  532.         if ((P[Y] < Prism->Height1) || (P[Y] > Prism->Height2))
  533.         {
  534.           return(FALSE);
  535.         }
  536.       }
  537.       else
  538.       {
  539.         /* Intersect ray with the cap-plane. */
  540.  
  541.         k = (Prism->Height2 - P[Y]) / D[Y];
  542.  
  543.         if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  544.         {
  545.           u = (P[X] + k * D[X]) / Prism->Height2;
  546.           v = (P[Z] + k * D[Z]) / Prism->Height2;
  547.  
  548.           if (in_curve(Prism, u, v))
  549.           {
  550.             Intersection[i].t   = CAP_HIT;
  551.             Intersection[i++].d = k / len;
  552.           }
  553.         }
  554.  
  555.         /* Intersect ray with the base-plane. */
  556.  
  557.         if (Prism->Height1 > 0.0)
  558.         {
  559.           k = (Prism->Height1 - P[Y]) / D[Y];
  560.  
  561.           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  562.           {
  563.             u = (P[X] + k * D[X]) / Prism->Height1;
  564.             v = (P[Z] + k * D[Z]) / Prism->Height1;
  565.  
  566.             if (in_curve(Prism, u, v))
  567.             {
  568.               Intersection[i].t   = BASE_HIT;
  569.               Intersection[i++].d = k / len;
  570.             }
  571.           }
  572.         }
  573.       }
  574.  
  575.       /* Precompute ray-only dependant constants. */
  576.  
  577.       k1 = P[Z] * D[Y] - P[Y] * D[Z];
  578.  
  579.       k2 = P[Y] * D[X] - P[X] * D[Y];
  580.  
  581.       k3 = P[X] * D[Z] - P[Z] * D[X];
  582.  
  583.       /* Intersect ray with the spline segments. */
  584.  
  585.       if ((fabs(D[X]) > EPSILON) || (fabs(D[Z]) > EPSILON))
  586.       {
  587.         for (j = 0; j < Prism->Number; j++)
  588.         {
  589.           Entry = Prism->Spline->Entry[j];
  590.  
  591.           /* Test spline's bounding rectangle (modified Cohen-Sutherland). */
  592.  
  593.           if (((D[X] >= 0.0) && (P[X] > Entry.x2)) ||
  594.               ((D[X] <= 0.0) && (P[X] < Entry.x1)) ||
  595.               ((D[Z] >= 0.0) && (P[Z] > Entry.y2)) ||
  596.               ((D[Z] <= 0.0) && (P[Z] < Entry.y1)))
  597.           {
  598.             continue;
  599.           }
  600.  
  601.           /* Number of roots found. */
  602.  
  603.           n = 0;
  604.  
  605.           switch (Prism->Spline_Type)
  606.           {
  607.             case LINEAR_SPLINE :
  608.  
  609.               /* Solve linear equation. */
  610.  
  611.               x[0] = Entry.C[X] * k1 + Entry.C[Y] * k2;
  612.  
  613.               x[1] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
  614.  
  615.               if (fabs(x[0]) > EPSILON)
  616.               {
  617.                 y[n++] = -x[1] / x[0];
  618.               }
  619.  
  620.               break;
  621.  
  622.             case QUADRATIC_SPLINE :
  623.  
  624.               /* Solve quadratic equation. */
  625.  
  626.               x[0] = Entry.B[X] * k1 + Entry.B[Y] * k2;
  627.  
  628.               x[1] = Entry.C[X] * k1 + Entry.C[Y] * k2;
  629.  
  630.               x[2] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
  631.  
  632.               n = Solve_Polynomial(2, x, y, FALSE, 0.0);
  633.  
  634.               break;
  635.  
  636.             case CUBIC_SPLINE :
  637.  
  638.               /* Solve cubic equation. */
  639.  
  640.               x[0] = Entry.A[X] * k1 + Entry.A[Y] * k2;
  641.  
  642.               x[1] = Entry.B[X] * k1 + Entry.B[Y] * k2;
  643.  
  644.               x[2] = Entry.C[X] * k1 + Entry.C[Y] * k2;
  645.  
  646.               x[3] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
  647.  
  648.               n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
  649.  
  650.               break;
  651.           }
  652.  
  653.           /* Test roots for valid intersections. */
  654.  
  655.           while (n--)
  656.           {
  657.             w = y[n];
  658.  
  659.             if ((w >= 0.0) && (w <= 1.0))
  660.             {
  661.               k = w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X];
  662.  
  663.               h = D[X] - k * D[Y];
  664.  
  665.               if (fabs(h) > EPSILON)
  666.               {
  667.                 k = (k * P[Y] - P[X]) / h;
  668.               }
  669.               else
  670.               {
  671.                 k = w * (w * (w * Entry.A[Y] + Entry.B[Y]) + Entry.C[Y]) + Entry.D[Y];
  672.  
  673.                 h = D[Z] - k * D[Y];
  674.  
  675.                 if (fabs(h) > EPSILON)
  676.                 {
  677.                   k = (k * P[Y] - P[Z]) / h;
  678.                 }
  679.                 else
  680.                 {
  681.                   /* This should never happen! */
  682.                   continue;
  683.                 }
  684.               }
  685.  
  686.               /* Verify that intersection height is valid. */
  687.  
  688.               h = P[Y] + k * D[Y];
  689.  
  690.               if ((h >= Prism->Height1) && (h <= Prism->Height2))
  691.               {
  692.                 Intersection[i].t   = SPLINE_HIT;
  693.                 Intersection[i].n   = j;
  694.                 Intersection[i].w   = w;
  695.                 Intersection[i++].d = k / len;
  696.               }
  697.             }
  698.           }
  699.         }
  700.       }
  701.  
  702.       break;
  703.  
  704.       default:
  705.  
  706.         Error("Unknown sweep type in intersect_prism().\n");
  707.   }
  708.  
  709.   if (i)
  710.   {
  711.     Increase_Counter(stats[Ray_Prism_Tests_Succeeded]);
  712.   }
  713.  
  714.   return(i);
  715. }
  716.  
  717.  
  718.  
  719. /*****************************************************************************
  720. *
  721. * FUNCTION
  722. *
  723. *   Inside_Prism
  724. *
  725. * INPUT
  726. *
  727. *   IPoint - Intersection point
  728. *   Object - Object
  729. *   
  730. * OUTPUT
  731. *   
  732. * RETURNS
  733. *
  734. *   int - TRUE if inside
  735. *   
  736. * AUTHOR
  737. *
  738. *   Dieter Bayer
  739. *   
  740. * DESCRIPTION
  741. *
  742. *   Test if point lies inside a prism.
  743. *
  744. * CHANGES
  745. *
  746. *   May 1994 : Creation.
  747. *
  748. ******************************************************************************/
  749.  
  750. static int Inside_Prism(IPoint, Object)
  751. VECTOR IPoint;
  752. OBJECT *Object;
  753. {
  754.   VECTOR P;
  755.   PRISM *Prism = (PRISM *)Object;
  756.  
  757.   /* Transform the point into the prism space. */
  758.  
  759.   MInvTransPoint(P, IPoint, Prism->Trans);
  760.  
  761.   if ((P[Y] >= Prism->Height1) && (P[Y] < Prism->Height2))
  762.   {
  763.     if (Prism->Sweep_Type == CONIC_SWEEP)
  764.     {
  765.       /* Scale x and z coordinate. */
  766.  
  767.       if (P[Y] > 0.0)
  768.       {
  769.         P[X] /= P[Y];
  770.         P[Z] /= P[Y];
  771.       }
  772.       else
  773.       {
  774.         P[X] = P[Z] = HUGE_VAL;
  775.       }
  776.     }
  777.  
  778.     if (in_curve(Prism, P[X], P[Z]))
  779.     {
  780.       return(!Test_Flag(Prism, INVERTED_FLAG));
  781.     }
  782.   }
  783.  
  784.   return(Test_Flag(Prism, INVERTED_FLAG));
  785. }
  786.  
  787.  
  788.  
  789. /*****************************************************************************
  790. *
  791. * FUNCTION
  792. *
  793. *   Prism_Normal
  794. *
  795. * INPUT
  796. *
  797. *   Result - Normal vector
  798. *   Object - Object
  799. *   Inter  - Intersection found
  800. *   
  801. * OUTPUT
  802. *
  803. *   Result
  804. *   
  805. * RETURNS
  806. *   
  807. * AUTHOR
  808. *
  809. *   Dieter Bayer
  810. *   
  811. * DESCRIPTION
  812. *
  813. *   Calculate the normal of the prism in a given point.
  814. *
  815. * CHANGES
  816. *
  817. *   May 1994 : Creation.
  818. *
  819. ******************************************************************************/
  820.  
  821. static void Prism_Normal(Result, Object, Inter)
  822. OBJECT *Object;
  823. VECTOR Result;
  824. INTERSECTION *Inter;
  825. {
  826.   DBL r;
  827.   VECTOR P;
  828.   PRISM_SPLINE_ENTRY Entry;
  829.   PRISM *Prism = (PRISM *)Object;
  830.   VECTOR N;
  831.  
  832.   Make_Vector(N, 0.0, 1.0, 0.0);
  833.  
  834.   if (Inter->i1 == SPLINE_HIT)
  835.   {
  836.     Entry = Prism->Spline->Entry[Inter->i2];
  837.  
  838.     switch (Prism->Sweep_Type)
  839.     {
  840.       case LINEAR_SWEEP:
  841.  
  842.         N[X] =   Inter->d1 * (3.0 * Entry.A[Y] * Inter->d1 + 2.0 * Entry.B[Y]) + Entry.C[Y];
  843.         N[Y] =   0.0;
  844.         N[Z] = -(Inter->d1 * (3.0 * Entry.A[X] * Inter->d1 + 2.0 * Entry.B[X]) + Entry.C[X]);
  845.  
  846.         break;
  847.  
  848.       case CONIC_SWEEP:
  849.  
  850.         /* Transform the point into the prism space. */
  851.  
  852.         MInvTransPoint(P, Inter->IPoint, Prism->Trans);
  853.  
  854.         if (P[Y] > 0.0)
  855.         {
  856.           r = sqrt(P[X] * P[X] + P[Z] * P[Z]);
  857.  
  858.           N[X] =  Inter->d1 * (3.0 * Entry.A[Y] * Inter->d1 + 2.0 * Entry.B[Y]) + Entry.C[Y];
  859.           N[Y] = -r / P[Y];
  860.           N[Z] = -(Inter->d1 * (3.0 * Entry.A[X] * Inter->d1 + 2.0 * Entry.B[X]) + Entry.C[X]);
  861.         }
  862.  
  863.         break;
  864.  
  865.       default:
  866.  
  867.         Error("Unknown sweep type in Prism_Normal().\n");
  868.     }
  869.   }
  870.  
  871.   /* Transform the normalt out of the prism space. */
  872.  
  873.   MTransNormal(Result, N, Prism->Trans);
  874.  
  875.   VNormalize(Result, Result);
  876. }
  877.  
  878.  
  879.  
  880. /*****************************************************************************
  881. *
  882. * FUNCTION
  883. *
  884. *   Translate_Prism
  885. *
  886. * INPUT
  887. *
  888. *   Object - Object
  889. *   Vector - Translation vector
  890. *   
  891. * OUTPUT
  892. *
  893. *   Object
  894. *   
  895. * RETURNS
  896. *   
  897. * AUTHOR
  898. *
  899. *   Dieter Bayer
  900. *   
  901. * DESCRIPTION
  902. *
  903. *   Translate a prism.
  904. *
  905. * CHANGES
  906. *
  907. *   May 1994 : Creation.
  908. *
  909. ******************************************************************************/
  910.  
  911. static void Translate_Prism(Object, Vector, Trans)
  912. OBJECT *Object;
  913. VECTOR Vector;
  914. TRANSFORM *Trans;
  915. {
  916.   Transform_Prism(Object, Trans);
  917. }
  918.  
  919.  
  920.  
  921. /*****************************************************************************
  922. *
  923. * FUNCTION
  924. *
  925. *   Rotate_Prism
  926. *
  927. * INPUT
  928. *
  929. *   Object - Object
  930. *   Vector - Rotation vector
  931. *   
  932. * OUTPUT
  933. *
  934. *   Object
  935. *   
  936. * RETURNS
  937. *   
  938. * AUTHOR
  939. *
  940. *   Dieter Bayer
  941. *   
  942. * DESCRIPTION
  943. *
  944. *   Rotate a prism.
  945. *
  946. * CHANGES
  947. *
  948. *   May 1994 : Creation.
  949. *
  950. ******************************************************************************/
  951.  
  952. static void Rotate_Prism(Object, Vector, Trans)
  953. OBJECT *Object;
  954. VECTOR Vector;
  955. TRANSFORM *Trans;
  956. {
  957.   Transform_Prism(Object, Trans);
  958. }
  959.  
  960.  
  961.  
  962. /*****************************************************************************
  963. *
  964. * FUNCTION
  965. *
  966. *   Scale_Prism
  967. *
  968. * INPUT
  969. *
  970. *   Object - Object
  971. *   Vector - Scaling vector
  972. *   
  973. * OUTPUT
  974. *
  975. *   Object
  976. *   
  977. * RETURNS
  978. *   
  979. * AUTHOR
  980. *
  981. *   Dieter Bayer
  982. *
  983. * DESCRIPTION
  984. *
  985. *   Scale a prism.
  986. *
  987. * CHANGES
  988. *
  989. *   May 1994 : Creation.
  990. *
  991. ******************************************************************************/
  992.  
  993. static void Scale_Prism(Object, Vector, Trans)
  994. OBJECT *Object;
  995. VECTOR Vector;
  996. TRANSFORM *Trans;
  997. {
  998.   Transform_Prism(Object, Trans);
  999. }
  1000.  
  1001.  
  1002.  
  1003. /*****************************************************************************
  1004. *
  1005. * FUNCTION
  1006. *
  1007. *   Transform_Prism
  1008. *
  1009. * INPUT
  1010. *
  1011. *   Object - Object
  1012. *   Trans  - Transformation to apply
  1013. *   
  1014. * OUTPUT
  1015. *
  1016. *   Object
  1017. *   
  1018. * RETURNS
  1019. *   
  1020. * AUTHOR
  1021. *
  1022. *   Dieter Bayer
  1023. *   
  1024. * DESCRIPTION
  1025. *
  1026. *   Transform a prism and recalculate its bounding box.
  1027. *
  1028. * CHANGES
  1029. *
  1030. *   May 1994 : Creation.
  1031. *
  1032. ******************************************************************************/
  1033.  
  1034. static void Transform_Prism(Object, Trans)
  1035. OBJECT *Object;
  1036. TRANSFORM *Trans;
  1037. {
  1038.   Compose_Transforms(((PRISM *)Object)->Trans, Trans);
  1039.  
  1040.   Compute_Prism_BBox((PRISM *)Object);
  1041. }
  1042.  
  1043.  
  1044.  
  1045. /*****************************************************************************
  1046. *
  1047. * FUNCTION
  1048. *
  1049. *   Invert_Prism
  1050. *
  1051. * INPUT
  1052. *
  1053. *   Object - Object
  1054. *   
  1055. * OUTPUT
  1056. *
  1057. *   Object
  1058. *   
  1059. * RETURNS
  1060. *   
  1061. * AUTHOR
  1062. *
  1063. *   Dieter Bayer
  1064. *   
  1065. * DESCRIPTION
  1066. *
  1067. *   Invert a prism.
  1068. *
  1069. * CHANGES
  1070. *
  1071. *   May 1994 : Creation.
  1072. *
  1073. ******************************************************************************/
  1074.  
  1075. static void Invert_Prism(Object)
  1076. OBJECT *Object;
  1077. {
  1078.   Invert_Flag(Object, INVERTED_FLAG);
  1079. }
  1080.  
  1081.  
  1082.  
  1083. /*****************************************************************************
  1084. *
  1085. * FUNCTION
  1086. *
  1087. *   Create_Prism
  1088. *
  1089. * INPUT
  1090. *   
  1091. * OUTPUT
  1092. *   
  1093. * RETURNS
  1094. *
  1095. *   PRISM * - new prism
  1096. *   
  1097. * AUTHOR
  1098. *
  1099. *   Dieter Bayer
  1100. *   
  1101. * DESCRIPTION
  1102. *
  1103. *   Create a new prism.
  1104. *
  1105. * CHANGES
  1106. *
  1107. *   May 1994 : Creation.
  1108. *
  1109. ******************************************************************************/
  1110.  
  1111. PRISM *Create_Prism()
  1112. {
  1113.   PRISM *New;
  1114.  
  1115.   New = (PRISM *)POV_MALLOC(sizeof(PRISM), "prism");
  1116.  
  1117.   INIT_OBJECT_FIELDS(New,PRISM_OBJECT,&Prism_Methods)
  1118.  
  1119.   New->Trans = Create_Transform();
  1120.  
  1121.   New->x1      =
  1122.   New->x2      =
  1123.   New->y1      =
  1124.   New->y2      =
  1125.   New->Height1 =
  1126.   New->Height2 = 0.0;
  1127.  
  1128.   New->Number = 0;
  1129.  
  1130.   New->Spline_Type = LINEAR_SPLINE;
  1131.   New->Sweep_Type  = LINEAR_SWEEP;
  1132.  
  1133.   New->Spline = NULL;
  1134.  
  1135.   New->Intersections = NULL;
  1136.  
  1137.   Set_Flag(New, CLOSED_FLAG);
  1138.  
  1139.   return(New);
  1140. }
  1141.  
  1142.  
  1143.  
  1144. /*****************************************************************************
  1145. *
  1146. * FUNCTION
  1147. *
  1148. *   Copy_Prism
  1149. *
  1150. * INPUT
  1151. *
  1152. *   Object - Object
  1153. *   
  1154. * OUTPUT
  1155. *   
  1156. * RETURNS
  1157. *
  1158. *   void * - New prism
  1159. *   
  1160. * AUTHOR
  1161. *
  1162. *   Dieter Bayer
  1163. *   
  1164. * DESCRIPTION
  1165. *
  1166. *   Copy a prism structure.
  1167. *
  1168. *   NOTE: The splines are not copied, only the number of references is
  1169. *         counted, so that Destray_Prism() knows if they can be destroyed.
  1170. *
  1171. * CHANGES
  1172. *
  1173. *   May 1994 : Creation.
  1174. *
  1175. *   Sep 1994 : fixed memory leakage [DB]
  1176. *
  1177. ******************************************************************************/
  1178.  
  1179. static void *Copy_Prism(Object)
  1180. OBJECT *Object;
  1181. {
  1182.   PRISM *New, *Prism = (PRISM *)Object;
  1183.  
  1184.   New = Create_Prism();
  1185.  
  1186.   /* Get rid of the transformation created in Create_Prism(). */
  1187.  
  1188.   Destroy_Transform(New->Trans);
  1189.  
  1190.   /* Copy prism. */
  1191.  
  1192.   *New = *Prism;
  1193.  
  1194.   New->Trans = Copy_Transform(((PRISM *)Object)->Trans);
  1195.  
  1196.   New->Spline->References++;
  1197.  
  1198.   Prism->Intersections = (PRISM_INT *)POV_MALLOC((Prism->Number+2)*sizeof(PRISM_INT), "prism intersection list");
  1199.  
  1200.   return(New);
  1201. }
  1202.  
  1203.  
  1204.  
  1205. /*****************************************************************************
  1206. *
  1207. * FUNCTION
  1208. *
  1209. *   Destroy_Prism
  1210. *
  1211. * INPUT
  1212. *
  1213. *   Object - Object
  1214. *   
  1215. * OUTPUT
  1216. *
  1217. *   Object
  1218. *   
  1219. * RETURNS
  1220. *   
  1221. * AUTHOR
  1222. *
  1223. *   Dieter Bayer
  1224. *   
  1225. * DESCRIPTION
  1226. *
  1227. *   Destroy a prism.
  1228. *
  1229. *   NOTE: The splines are destroyed if they are no longer used by any copy.
  1230. *
  1231. * CHANGES
  1232. *
  1233. *   May 1994 : Creation.
  1234. *
  1235. ******************************************************************************/
  1236.  
  1237. static void Destroy_Prism (Object)
  1238. OBJECT *Object;
  1239. {
  1240.   PRISM *Prism = (PRISM *)Object;
  1241.  
  1242.   Destroy_Transform(Prism->Trans);
  1243.  
  1244.   POV_FREE(Prism->Intersections);
  1245.  
  1246.   if (--(Prism->Spline->References) == 0)
  1247.   {
  1248.     POV_FREE(Prism->Spline->Entry);
  1249.  
  1250.     POV_FREE(Prism->Spline);
  1251.   }
  1252.  
  1253.   POV_FREE(Object);
  1254. }
  1255.  
  1256.  
  1257.  
  1258. /*****************************************************************************
  1259. *
  1260. * FUNCTION
  1261. *
  1262. *   Compute_Prism_BBox
  1263. *
  1264. * INPUT
  1265. *
  1266. *   Prism - Prism
  1267. *   
  1268. * OUTPUT
  1269. *
  1270. *   Prism
  1271. *   
  1272. * RETURNS
  1273. *   
  1274. * AUTHOR
  1275. *
  1276. *   Dieter Bayer
  1277. *   
  1278. * DESCRIPTION
  1279. *
  1280. *   Calculate the bounding box of a prism.
  1281. *
  1282. * CHANGES
  1283. *
  1284. *   May 1994 : Creation.
  1285. *
  1286. ******************************************************************************/
  1287.  
  1288. void Compute_Prism_BBox(Prism)
  1289. PRISM *Prism;
  1290. {
  1291.   Make_BBox(Prism->BBox, Prism->x1, Prism->Height1, Prism->y1,
  1292.     Prism->x2 - Prism->x1, Prism->Height2 - Prism->Height1, Prism->y2 - Prism->y1);
  1293.  
  1294.   Recompute_BBox(&Prism->BBox, Prism->Trans);
  1295. }
  1296.  
  1297.  
  1298.  
  1299. /*****************************************************************************
  1300. *
  1301. * FUNCTION
  1302. *
  1303. *   in_curve
  1304. *
  1305. * INPUT
  1306. *
  1307. *   Prism - Prism to test
  1308. *   u, v  - Coordinates
  1309. *   
  1310. * OUTPUT
  1311. *   
  1312. * RETURNS
  1313. *
  1314. *   int - TRUE if inside
  1315. *   
  1316. * AUTHOR
  1317. *
  1318. *   Dieter Bayer, June 1994
  1319. *   
  1320. * DESCRIPTION
  1321. *
  1322. *   Test if a 2d point lies inside a prism's spline curve.
  1323. *
  1324. *   We travel from the current point in positive u direction
  1325. *   and count the number of crossings with the curve.
  1326. *
  1327. * CHANGES
  1328. *
  1329. *   May 1994 : Creation.
  1330. *
  1331. ******************************************************************************/
  1332.  
  1333. static int in_curve(Prism, u, v)
  1334. PRISM *Prism;
  1335. DBL u, v;
  1336. {
  1337.   int i, n, NC;
  1338.   DBL k, w;
  1339.   DBL x[4], y[3];
  1340.   PRISM_SPLINE_ENTRY Entry;
  1341.  
  1342.   NC = 0;
  1343.  
  1344.   /* First test overall bounding rectangle. */
  1345.   
  1346.   if ((u >= Prism->x1) && (u <= Prism->x2) &&
  1347.       (v >= Prism->y1) && (v <= Prism->y2))
  1348.   {
  1349.     for (i = 0; i < Prism->Number; i++)
  1350.     {
  1351.       Entry = Prism->Spline->Entry[i];
  1352.  
  1353.       /* Test if current segment can be hit. */
  1354.  
  1355.       if ((v >= Entry.y1) && (v <= Entry.y2) && (u <= Entry.x2))
  1356.       {
  1357.         x[0] = Entry.A[Y];
  1358.         x[1] = Entry.B[Y];
  1359.         x[2] = Entry.C[Y];
  1360.         x[3] = Entry.D[Y] - v;
  1361.  
  1362.         n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
  1363.  
  1364.         while (n--)
  1365.         {
  1366.           w = y[n];
  1367.  
  1368.           if ((w >= 0.0) && (w <= 1.0))
  1369.           {
  1370.             k  = w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X] - u;
  1371.  
  1372.             if (k >= 0.0)
  1373.             {
  1374.               NC++;
  1375.             }
  1376.           }
  1377.         }
  1378.       }
  1379.     }
  1380.   }
  1381.  
  1382.   return(NC & 1);
  1383. }
  1384.  
  1385.  
  1386.  
  1387. /*****************************************************************************
  1388. *
  1389. * FUNCTION
  1390. *
  1391. *   test_rectangle
  1392. *
  1393. * INPUT
  1394. *   
  1395. * OUTPUT
  1396. *   
  1397. * RETURNS
  1398. *   
  1399. * AUTHOR
  1400. *
  1401. *   Dieter Bayer, July 1994
  1402. *   
  1403. * DESCRIPTION
  1404. *
  1405. *   Test if the 2d ray (= P + t * D) intersects a rectangle.
  1406. *
  1407. * CHANGES
  1408. *
  1409. *   May 1994 : Creation.
  1410. *
  1411. ******************************************************************************/
  1412.  
  1413. static int test_rectangle(P, D, x1, z1, x2, z2)
  1414. VECTOR P, D;
  1415. DBL x1, z1, x2, z2;
  1416. {
  1417.   DBL dmin, dmax, tmin, tmax;
  1418.  
  1419.   if (fabs(D[X]) > EPSILON)
  1420.   {
  1421.     if (D[X] > 0.0)
  1422.     {
  1423.       dmin = (x1 - P[X]) / D[X];
  1424.       dmax = (x2 - P[X]) / D[X];
  1425.  
  1426.       if (dmax < EPSILON)
  1427.       {
  1428.         return(FALSE);
  1429.       }
  1430.     }
  1431.     else
  1432.     {
  1433.       dmax = (x1 - P[X]) / D[X];
  1434.  
  1435.       if (dmax < EPSILON)
  1436.       {
  1437.         return(FALSE);
  1438.       }
  1439.  
  1440.       dmin = (x2 - P[X]) / D[X];
  1441.     }
  1442.  
  1443.     if (dmin > dmax)
  1444.     {
  1445.       return(FALSE);
  1446.     }
  1447.   }
  1448.   else
  1449.   {
  1450.     if ((P[X] < x1) || (P[X] > x2))
  1451.     {
  1452.       return(FALSE);
  1453.     }
  1454.     else
  1455.     {
  1456.       dmin = -BOUND_HUGE;
  1457.       dmax =  BOUND_HUGE;
  1458.     }
  1459.   }
  1460.  
  1461.   if (fabs(D[Z]) > EPSILON)
  1462.   {
  1463.     if (D[Z] > 0.0)
  1464.     {
  1465.       tmin = (z1 - P[Z]) / D[Z];
  1466.       tmax = (z2 - P[Z]) / D[Z];
  1467.     }
  1468.     else
  1469.     {
  1470.       tmax = (z1 - P[Z]) / D[Z];
  1471.       tmin = (z2 - P[Z]) / D[Z];
  1472.     }
  1473.  
  1474.     if (tmax < dmax)
  1475.     {
  1476.       if (tmax < EPSILON)
  1477.       {
  1478.         return(FALSE);
  1479.       }
  1480.  
  1481.       if (tmin > dmin)
  1482.       {
  1483.         if (tmin > tmax)
  1484.         {
  1485.           return(FALSE);
  1486.         }
  1487.  
  1488.         dmin = tmin;
  1489.       }
  1490.       else
  1491.       {
  1492.         if (dmin > tmax)
  1493.         {
  1494.           return(FALSE);
  1495.         }
  1496.       }
  1497.  
  1498.       /*dmax = tmax; */ /*not needed CEY[1/95]*/
  1499.     }
  1500.     else
  1501.     {
  1502.       if (tmin > dmin)
  1503.       {
  1504.         if (tmin > dmax)
  1505.         {
  1506.           return(FALSE);
  1507.         }
  1508.  
  1509.         /* dmin = tmin; */  /*not needed CEY[1/95]*/
  1510.       }
  1511.     }
  1512.   }
  1513.   else
  1514.   {
  1515.     if ((P[Z] < z1) || (P[Z] > z2))
  1516.     {
  1517.       return(FALSE);
  1518.     }
  1519.   }
  1520.  
  1521.   return(TRUE);
  1522. }
  1523.  
  1524.  
  1525.  
  1526. /*****************************************************************************
  1527. *
  1528. * FUNCTION
  1529. *
  1530. *   Compute_Prism
  1531. *
  1532. * INPUT
  1533. *
  1534. *   Prism - Prism
  1535. *   P     - Points defining prism
  1536. *   
  1537. * OUTPUT
  1538. *
  1539. *   Prism
  1540. *
  1541. * RETURNS
  1542. *
  1543. * AUTHOR
  1544. *
  1545. *   Dieter Bayer, June 1994
  1546. *
  1547. * DESCRIPTION
  1548. *
  1549. *   Calculate the spline segments of a prism from a set of points.
  1550. *
  1551. * CHANGES
  1552. *
  1553. *   May 1994 : Creation.
  1554. *
  1555. ******************************************************************************/
  1556.  
  1557. void Compute_Prism(Prism, P)
  1558. PRISM *Prism;
  1559. UV_VECT *P;
  1560. {
  1561.   int i, n, number_of_splines;
  1562.   int i1, i2, i3;
  1563.  
  1564.   DBL x[4], xmin, xmax;
  1565.   DBL y[4], ymin, ymax;
  1566.   DBL c[3], r[2];
  1567.  
  1568.   UV_VECT A, B, C, D, First;
  1569.  
  1570.   /* Allocate Object->Number segments. */
  1571.  
  1572.   if (Prism->Spline == NULL)
  1573.   {
  1574.     Prism->Spline = (PRISM_SPLINE *)POV_MALLOC(sizeof(PRISM_SPLINE), "spline segments of prism");
  1575.  
  1576.     Prism->Spline->References = 1;
  1577.  
  1578.     Prism->Spline->Entry = (PRISM_SPLINE_ENTRY *)POV_MALLOC(Prism->Number*sizeof(PRISM_SPLINE_ENTRY), "spline segments of prism");
  1579.   }
  1580.   else
  1581.   {
  1582.     /* This should never happen! */
  1583.  
  1584.     Error("Prism segments are already defined.\n");
  1585.   }
  1586.  
  1587.   /* Allocate prism intersections list. */
  1588.  
  1589.   Prism->Intersections = (PRISM_INT *)POV_MALLOC((Prism->Number+2)*sizeof(PRISM_INT), "prism intersection list");
  1590.  
  1591.   /***************************************************************************
  1592.   * Calculate segments.
  1593.   ****************************************************************************/
  1594.  
  1595.   /* We want to know the size of the overall bounding rectangle. */
  1596.  
  1597.   xmax = ymax = -BOUND_HUGE;
  1598.   xmin = ymin =  BOUND_HUGE;
  1599.  
  1600.   /* Get first segment point. */
  1601.  
  1602.   switch (Prism->Spline_Type)
  1603.   {
  1604.     case LINEAR_SPLINE:
  1605.       Assign_UV_Vect(First, P[0]);
  1606.  
  1607.       break;
  1608.  
  1609.     case QUADRATIC_SPLINE:
  1610.     case CUBIC_SPLINE:
  1611.  
  1612.       Assign_UV_Vect(First, P[1]);
  1613.  
  1614.       break;
  1615.   }
  1616.  
  1617.   for (i = number_of_splines = 0; i < Prism->Number-1; i++)
  1618.   {
  1619.     /* Get indices of previous and next two points. */
  1620.  
  1621.     i1 = i + 1;
  1622.     i2 = i + 2;
  1623.     i3 = i + 3;
  1624.  
  1625.     switch (Prism->Spline_Type)
  1626.     {
  1627.       /*************************************************************************
  1628.       * Linear spline (nothing more than a simple polygon).
  1629.       **************************************************************************/
  1630.  
  1631.       case LINEAR_SPLINE :
  1632.  
  1633.         if (i1 >= Prism->Number)
  1634.         {
  1635.           Error("Too few points in prism. Prism not closed? Control points missing?\n");
  1636.         }
  1637.  
  1638.         /* Use linear interpolation. */
  1639.  
  1640.         A[X] =  0.0;
  1641.         B[X] =  0.0;
  1642.         C[X] = -1.0 * P[i][X] + 1.0 * P[i1][X];
  1643.         D[X] =  1.0 * P[i][X];
  1644.  
  1645.         A[Y] =  0.0;
  1646.         B[Y] =  0.0;
  1647.         C[Y] = -1.0 * P[i][Y] + 1.0 * P[i1][Y];
  1648.         D[Y] =  1.0 * P[i][Y];
  1649.  
  1650.         x[0] = x[2] = P[i][X];
  1651.         x[1] = x[3] = P[i1][X];
  1652.  
  1653.         y[0] = y[2] = P[i][Y];
  1654.         y[1] = y[3] = P[i1][Y];
  1655.  
  1656.         break;
  1657.  
  1658.       /*************************************************************************
  1659.       * Quadratic spline.
  1660.       **************************************************************************/
  1661.  
  1662.       case QUADRATIC_SPLINE :
  1663.  
  1664.         if (i2 >= Prism->Number)
  1665.         {
  1666.           Error("Too few points in prism.\n");
  1667.         }
  1668.  
  1669.         /* Use quadratic interpolation. */
  1670.  
  1671.         A[X] =  0.0;
  1672.         B[X] =  0.5 * P[i][X] - 1.0 * P[i1][X] + 0.5 * P[i2][X];
  1673.         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
  1674.         D[X] =                  1.0 * P[i1][X];
  1675.  
  1676.         A[Y] =  0.0;
  1677.         B[Y] =  0.5 * P[i][Y] - 1.0 * P[i1][Y] + 0.5 * P[i2][Y];
  1678.         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
  1679.         D[Y] =                  1.0 * P[i1][Y];
  1680.  
  1681.         x[0] = x[2] = P[i1][X];
  1682.         x[1] = x[3] = P[i2][X];
  1683.  
  1684.         y[0] = y[2] = P[i1][Y];
  1685.         y[1] = y[3] = P[i2][Y];
  1686.  
  1687.         break;
  1688.  
  1689.       /*************************************************************************
  1690.       * Cubic spline.
  1691.       **************************************************************************/
  1692.  
  1693.       case CUBIC_SPLINE :
  1694.  
  1695.         if (i3 >= Prism->Number)
  1696.         {
  1697.           Error("Too few points in prism.\n");
  1698.         }
  1699.  
  1700.         /* Use cubic interpolation. */
  1701.  
  1702.         A[X] = -0.5 * P[i][X] + 1.5 * P[i1][X] - 1.5 * P[i2][X] + 0.5 * P[i3][X];
  1703.         B[X] =        P[i][X] - 2.5 * P[i1][X] + 2.0 * P[i2][X] - 0.5 * P[i3][X];
  1704.         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
  1705.         D[X] =                        P[i1][X];
  1706.  
  1707.         A[Y] = -0.5 * P[i][Y] + 1.5 * P[i1][Y] - 1.5 * P[i2][Y] + 0.5 * P[i3][Y];
  1708.         B[Y] =        P[i][Y] - 2.5 * P[i1][Y] + 2.0 * P[i2][Y] - 0.5 * P[i3][Y];
  1709.         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
  1710.         D[Y] =                        P[i1][Y];
  1711.  
  1712.         x[0] = x[2] = P[i1][X];
  1713.         x[1] = x[3] = P[i2][X];
  1714.  
  1715.         y[0] = y[2] = P[i1][Y];
  1716.         y[1] = y[3] = P[i2][Y];
  1717.  
  1718.         break;
  1719.  
  1720.       default:
  1721.  
  1722.         Error("Unknown spline type in Compute_Prism().\n");
  1723.     }
  1724.  
  1725.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].A, A);
  1726.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].B, B);
  1727.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].C, C);
  1728.     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].D, D);
  1729.  
  1730.     /* Get maximum coordinates in current segment. */
  1731.  
  1732.     c[0] = 3.0 * A[X];
  1733.     c[1] = 2.0 * B[X];
  1734.     c[2] = C[X];
  1735.  
  1736.     n = Solve_Polynomial(2, c, r, FALSE, 0.0);
  1737.  
  1738.     while (n--)
  1739.     {
  1740.       if ((r[n] >= 0.0) && (r[n] <= 1.0))
  1741.       {
  1742.         x[n] = r[n] * (r[n] * (r[n] * A[X] + B[X]) + C[X]) + D[X];
  1743.       }
  1744.     }
  1745.  
  1746.     c[0] = 3.0 * A[Y];
  1747.     c[1] = 2.0 * B[Y];
  1748.     c[2] = C[Y];
  1749.  
  1750.     n = Solve_Polynomial(2, c, r, FALSE, 0.0);
  1751.  
  1752.     while (n--)
  1753.     {
  1754.       if ((r[n] >= 0.0) && (r[n] <= 1.0))
  1755.       {
  1756.         y[n] = r[n] * (r[n] * (r[n] * A[Y] + B[Y]) + C[Y]) + D[Y];
  1757.       }
  1758.     }
  1759.  
  1760.     /* Set current segment's bounding rectangle. */
  1761.  
  1762.     Prism->Spline->Entry[number_of_splines].x1 = min(min(x[0], x[1]), min(x[2], x[3]));
  1763.     Prism->Spline->Entry[number_of_splines].x2 = max(max(x[0], x[1]), max(x[2], x[3]));
  1764.  
  1765.     Prism->Spline->Entry[number_of_splines].y1 = min(min(y[0], y[1]), min(y[2], y[3]));
  1766.     Prism->Spline->Entry[number_of_splines].y2 = max(max(y[0], y[1]), max(y[2], y[3]));
  1767.  
  1768.     /* Keep track of overall bounding rectangle. */
  1769.  
  1770.     xmin = min(xmin, Prism->Spline->Entry[number_of_splines].x1);
  1771.     xmax = max(xmax, Prism->Spline->Entry[number_of_splines].x2);
  1772.  
  1773.     ymin = min(ymin, Prism->Spline->Entry[number_of_splines].y1);
  1774.     ymax = max(ymax, Prism->Spline->Entry[number_of_splines].y2);
  1775.  
  1776.     number_of_splines++;
  1777.  
  1778.     /* Check if end of sub-prism is reached. */
  1779.  
  1780.     switch (Prism->Spline_Type)
  1781.     {
  1782.       case LINEAR_SPLINE:
  1783.  
  1784.         if ((fabs(P[i1][X] - First[X]) < EPSILON) &&
  1785.             (fabs(P[i1][Y] - First[Y]) < EPSILON))
  1786.         {
  1787.           i++;
  1788.  
  1789.           if (i+1 < Prism->Number)
  1790.           {
  1791.             Assign_UV_Vect(First, P[i+1]);
  1792.           }
  1793.         }
  1794.  
  1795.         break;
  1796.  
  1797.       case QUADRATIC_SPLINE:
  1798.  
  1799.         if ((fabs(P[i2][X] - First[X]) < EPSILON) &&
  1800.             (fabs(P[i2][Y] - First[Y]) < EPSILON))
  1801.         {
  1802.           i += 2;
  1803.  
  1804.           if (i+2 < Prism->Number)
  1805.           {
  1806.             Assign_UV_Vect(First, P[i+2]);
  1807.           }
  1808.         }
  1809.  
  1810.         break;
  1811.  
  1812.       case CUBIC_SPLINE:
  1813.  
  1814.         if ((fabs(P[i2][X] - First[X]) < EPSILON) &&
  1815.             (fabs(P[i2][Y] - First[Y]) < EPSILON))
  1816.         {
  1817.           i += 3;
  1818.  
  1819.           if (i+2 < Prism->Number)
  1820.           {
  1821.             Assign_UV_Vect(First, P[i+2]);
  1822.           }
  1823.         }
  1824.  
  1825.         break;
  1826.  
  1827.     }
  1828.   }
  1829.  
  1830.   Prism->Number = number_of_splines;
  1831.  
  1832.   /* Set overall bounding rectangle. */
  1833.  
  1834.   Prism->x1 = xmin;
  1835.   Prism->x2 = xmax;
  1836.  
  1837.   Prism->y1 = ymin;
  1838.   Prism->y2 = ymax;
  1839.  
  1840.   if (Prism->Sweep_Type == CONIC_SWEEP)
  1841.   {
  1842.     /* Recalculate bounding rectangles. */
  1843.  
  1844.     for (i = 0; i < Prism->Number; i++)
  1845.     {
  1846.       x[0] = Prism->Spline->Entry[i].x1 * Prism->Height1;
  1847.       x[1] = Prism->Spline->Entry[i].x1 * Prism->Height2;
  1848.       x[2] = Prism->Spline->Entry[i].x2 * Prism->Height1;
  1849.       x[3] = Prism->Spline->Entry[i].x2 * Prism->Height2;
  1850.  
  1851.       xmin = min(min(x[0], x[1]), min(x[2], x[3]));
  1852.       xmax = max(max(x[0], x[1]), max(x[2], x[3]));
  1853.  
  1854.       Prism->Spline->Entry[i].x1 = xmin;
  1855.       Prism->Spline->Entry[i].x2 = xmax;
  1856.  
  1857.       y[0] = Prism->Spline->Entry[i].y1 * Prism->Height1;
  1858.       y[1] = Prism->Spline->Entry[i].y1 * Prism->Height2;
  1859.       y[2] = Prism->Spline->Entry[i].y2 * Prism->Height1;
  1860.       y[3] = Prism->Spline->Entry[i].y2 * Prism->Height2;
  1861.  
  1862.       ymin = min(min(y[0], y[1]), min(y[2], y[3]));
  1863.       ymax = max(max(y[0], y[1]), max(y[2], y[3]));
  1864.  
  1865.       Prism->Spline->Entry[i].y1 = ymin;
  1866.       Prism->Spline->Entry[i].y2 = ymax;
  1867.     }
  1868.  
  1869.     /* Recalculate overall bounding rectangle. */
  1870.  
  1871.     x[0] = Prism->x1 * Prism->Height1;
  1872.     x[1] = Prism->x1 * Prism->Height2;
  1873.     x[2] = Prism->x2 * Prism->Height1;
  1874.     x[3] = Prism->x2 * Prism->Height2;
  1875.  
  1876.     xmin = min(min(x[0], x[1]), min(x[2], x[3]));
  1877.     xmax = max(max(x[0], x[1]), max(x[2], x[3]));
  1878.  
  1879.     Prism->x1 = xmin;
  1880.     Prism->x2 = xmax;
  1881.  
  1882.     y[0] = Prism->y1 * Prism->Height1;
  1883.     y[1] = Prism->y1 * Prism->Height2;
  1884.     y[2] = Prism->y2 * Prism->Height1;
  1885.     y[3] = Prism->y2 * Prism->Height2;
  1886.  
  1887.     ymin = min(min(y[0], y[1]), min(y[2], y[3]));
  1888.     ymax = max(max(y[0], y[1]), max(y[2], y[3]));
  1889.  
  1890.     Prism->y1 = ymin;
  1891.     Prism->y2 = ymax;
  1892.   }
  1893. }
  1894.  
  1895.  
  1896.  
  1897.